{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Datashader is a general-purpose tool for rasterizing (and re-rasterizing) data of many different types. To make it easier to apply this general functionality to the particular domain of geoscience, Datashader provides a few geospatial-specific utilities as well:\n",
    "\n",
    "* [Project points](#Project-points)\n",
    "* [Generate terrain](#Generate-terrain)\n",
    "* [Hillshade](#Hillshade)\n",
    "* [Slope](#Slope)\n",
    "* [Aspect](#Aspect)\n",
    "* [Bump map](#Bump-map)\n",
    "* [NDVI](#NDVI)\n",
    "* [Mean](#Mean)\n",
    "* [Proximity](#Proximity)\n",
    "* [Viewshed](#Viewshed)\n",
    "* [Zonal Statistics](#Zonal-Statistics)\n",
    "\n",
    "This notebook explains each of these topics in turn. See also [GeoViews](http://geoviews.org), which is designed to work with Datashader to provide a large range of additional geospatial functionality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Project points\n",
    "\n",
    "You can use [GeoViews](http://geoviews.org) or the underlying [pyproj/proj.4](https://pyproj4.github.io/pyproj) libraries to perform arbitrary projections to and from a large number of different coordinate reference systems. However, for the common case of wanting to view data with latitude and longitude coordinates on top of a Web Mercator tile source such as Google Maps or OpenStreetMap, Datashader provides a simple self-contained utility `lnglat_to_meters(longitude, latitude)` to project your data once, before visualization.  For instance, if you have a dataframe with some latitude and longitude points stretching from San Diego, California to Bangor, Maine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, pandas as pd\n",
    "from datashader.utils import lnglat_to_meters\n",
    "\n",
    "San_Diego = 32.715, -117.1625\n",
    "Bangor = 44.8, -68.8\n",
    "n = 20\n",
    "\n",
    "df = pd.DataFrame(dict(longitude = np.linspace(San_Diego[1], Bangor[1], n),\n",
    "                       latitude  = np.linspace(San_Diego[0], Bangor[0], n)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then you can create new columns (or overwrite old ones) with the projected points in meters from the origin (Web Mercator coordinates):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.loc[:, 'x'], df.loc[:, 'y'] = lnglat_to_meters(df.longitude,df.latitude)\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new x and y coordinates aren't very useful for humans to read, but they can now be overlaid directly onto web map sources, which are labeled with latitude and longitude appropriately by Bokeh (via HoloViews) but are actually in Web Mercator coordinates internally:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "from holoviews.operation.datashader import datashade, spread\n",
    "from holoviews.element import tiles\n",
    "hv.extension('bokeh')\n",
    "\n",
    "pts = spread(datashade(hv.Points(df, ['x', 'y']), cmap=\"white\", width=300, height=100), px=3)\n",
    "\n",
    "tiles.EsriImagery() * pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are using GeoViews, you can get the same effect by calling [gv.operation.project](http://geoviews.org/user_guide/Projections.html#Explicitly-projecting).  With GeoViews, you can also declare your object to be in lon,lat coordinates natively (`from cartopy import crs ; gv.Points(df, ['longitude', 'latitude'], crs=crs.PlateCarree())`) and let GeoViews then reproject the points as needed, but dynamic reprojection will be much slower for interactive use than projecting them in bulk ahead of time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"background-color: #f8d7da; border-radius: 5px; color: #721c24; padding: 10px\">Warning: As of datashader version 0.11.0 the datashader.geo module has been deprecated, and the functionality has migrated to <a href=\"https://github.com/makepath/xarray-spatial\">xarray-spatial</a>.</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Terrain Data\n",
    "\n",
    "The rest of the geo-related functions focus on raster data (or rasterized data, after a previous Datashader step that returns an Xarray object). To demonstrate using these raster-based functions, let's generate some fake terrain as an elevation raster:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, datashader as ds, datashader.geo as dsgeo\n",
    "from datashader.transfer_functions import shade, stack\n",
    "from datashader.colors import Elevation\n",
    "\n",
    "W = 800\n",
    "H = 600\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))\n",
    "terrain = dsgeo.generate_terrain(cvs)\n",
    "\n",
    "shade(terrain, cmap=['black', 'white'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grayscale value above shows the elevation linearly in intensity (with the large black areas indicating low elevation), but it will look more like a landscape if we map the lowest values to colors representing water, and the highest to colors representing mountaintops:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shade(terrain, cmap=Elevation, how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hillshade\n",
    "\n",
    "[Hillshade](https://en.wikipedia.org/wiki/Terrain_cartography) is a technique used to visualize terrain as shaded relief, illuminating it with a hypothetical light source. The illumination value for each cell is determined by its orientation to the light source, which is based on slope and aspect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "illuminated = dsgeo.hillshade(terrain)\n",
    "\n",
    "shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can combine hillshading with elevation colormapping to convey differences in terrain with elevation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack(shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear'),\n",
    "      shade(terrain,     cmap=Elevation,         alpha=128, how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slope\n",
    "[Slope](https://en.wikipedia.org/wiki/Slope) is the inclination of a surface. \n",
    "In geography, *slope* is amount of change in elevation of a terrain regarding its surroundings.\n",
    "\n",
    "Datashader's slope function returns slope in degrees.  Below we highlight areas at risk for avalanche by looking at [slopes around 38 degrees](http://wenatcheeoutdoors.org/2016/04/07/avalanche-abcs-for-snowshoers/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "risky = dsgeo.slope(terrain)\n",
    "risky.data = np.where(np.logical_and(risky.data > 25, risky.data < 50), 1, np.nan)\n",
    "\n",
    "stack(shade(terrain,     cmap=['black', 'white'], how='linear'),\n",
    "      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(risky,       cmap='red',              how='linear', alpha=200))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aspect\n",
    "\n",
    "[Aspect](https://en.wikipedia.org/wiki/Aspect_%28geography%29) is the orientation of slope, measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.\n",
    "\n",
    "Below, we look to find slopes that face close to North."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "north_faces = dsgeo.aspect(terrain)\n",
    "north_faces.data = np.where(np.logical_or(north_faces.data > 350 ,\n",
    "                                          north_faces.data < 10), 1, np.nan)\n",
    "\n",
    "stack(shade(terrain,     cmap=['black', 'white'], how='linear'),\n",
    "      shade(illuminated, cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(north_faces, cmap=['aqua'],           how='linear', alpha=100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NDVI\n",
    "\n",
    "The [Normalized Difference Vegetation Index](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) (NDVI) is a metric designed to detect regions with vegetation by measuring the difference between near-infrared (NIR) light (which vegetation reflects) and red light (which vegetation absorbs).\n",
    "\n",
    "The NDVI ranges over [-1,+1], where `-1` means more \"Red\" radiation while `+1` means more \"NIR\" radiation. NDVI values close to +1.0 suggest areas dense with active green foliage, while strongly negative values suggest cloud cover or snow, and values near zero suggest open water, urban areas, or bare soil. \n",
    "\n",
    "For our synthetic example here, we don't have access to NIR measurements, but we can approximate the results for demonstration purposes by using the green and blue channels of a colormapped image, as those represent a difference in wavelength similar to NIR vs. Red."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "\n",
    "rgba = stack(shade(terrain, cmap=Elevation, how='linear')).to_pil()\n",
    "r,g,b,a = [xr.DataArray(np.flipud(np.asarray(rgba.getchannel(c))))/255.0 \n",
    "           for c in ['R','G','B','A']]\n",
    "\n",
    "ndvi = dsgeo.ndvi(nir_agg=g, red_agg=b)\n",
    "shade(ndvi, cmap=['purple','black','green'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bump\n",
    "\n",
    "Bump mapping is a cartographic technique that can be used to create the appearance of trees or other land features, which is useful when synthesizing human-interpretable images from source data like land use classifications.\n",
    "\n",
    "`dsgeo.bump` will produce a bump aggregate for adding detail to the terrain.\n",
    "\n",
    "In this example, we will pretend the bumps are trees, and shade them with green.  We'll also use the elevation data to modulate whether there are trees and if so how tall they are.\n",
    "\n",
    "- First, we'll define a custom `height` function to return tree heights suitable for the given elevation range\n",
    "- `dsgeo.bump` accepts a function with only a single argument (`locations`), so we will use `functools.partial` to provide values for the other arguments.\n",
    "- Bump mapping isn't normally a performance bottleneck, but if you want, you can speed it up by using Numba on your custom `height` function (`from datashader.utils import ngjit`, then put `@ngjit` above `def heights(...)`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "def heights(locations, src, src_range, height=20):\n",
    "    num_bumps = locations.shape[0]\n",
    "    out = np.zeros(num_bumps, dtype=np.uint16)\n",
    "    for r in range(0, num_bumps):\n",
    "        loc = locations[r]\n",
    "        x = loc[0]\n",
    "        y = loc[1]\n",
    "        val = src[y, x]\n",
    "        if val >= src_range[0] and val < src_range[1]:\n",
    "            out[r] = height\n",
    "    return out\n",
    "\n",
    "T = 300000 # Number of trees to add per call\n",
    "src = terrain.data\n",
    "%time trees  = dsgeo.bump(W, H, count=T,    height_func=partial(heights, src=src, src_range=(1000, 1300), height=5))\n",
    "trees       += dsgeo.bump(W, H, count=T//2, height_func=partial(heights, src=src, src_range=(1300, 1700), height=20))\n",
    "trees       += dsgeo.bump(W, H, count=T//3, height_func=partial(heights, src=src, src_range=(1700, 2000), height=5))\n",
    "\n",
    "tree_colorize = trees.copy()\n",
    "tree_colorize.data[tree_colorize.data == 0] = np.nan\n",
    "hillshade = dsgeo.hillshade(terrain + trees)\n",
    "\n",
    "stack(shade(terrain,        cmap=['black', 'white'], how='linear'),\n",
    "      shade(hillshade,      cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(tree_colorize,  cmap='limegreen',        how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean\n",
    "The `datashader.mean` function will smooth a given aggregate by using a 3x3 mean convolution filter. Optional parameters include `passes`, which is used to run the mean filter multiple times, and also `excludes` which are values that will not be modified by the mean filter.\n",
    "\n",
    "We can use `mean` to add a coastal vignette to give out terrain scene a bit more character. Notice the water below now has a nice coastal gradient which adds some realism to our scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "LAND_CONSTANT = 50.0\n",
    "\n",
    "water = terrain.copy()\n",
    "water.data = np.where(water.data > 0, LAND_CONSTANT, 0)\n",
    "water = dsgeo.mean(water, passes=50, excludes=[LAND_CONSTANT])\n",
    "water.data[water.data == LAND_CONSTANT] = np.nan\n",
    "\n",
    "stack(shade(terrain,    cmap=['black', 'white'], how='linear'),\n",
    "      shade(water,      cmap=['aqua',  'white']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full scene\n",
    "\n",
    "We've now seen several of datashader's `geo` helper functions for working with elevation rasters.\n",
    "\n",
    "Let's make a full archipelago scene by stacking `terrain`, `water`, `hillshade`, and `tree_colorize` together into one output image: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack(shade(terrain,       cmap=Elevation,          how='linear'),\n",
    "      shade(water,         cmap=['aqua','white']),\n",
    "      shade(dsgeo.hillshade(terrain + trees),   cmap=['black', 'white'], how='linear', alpha=128),\n",
    "      shade(tree_colorize, cmap='limegreen',        how='linear'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Proximity\n",
    "\n",
    "The `datashader.spatial.proximity` function operates on a given aggregate to produce a new distance aggregate based on target values and a distance metric. The values in the new aggregate will be the distance (according to the given metric) between each array cell (pixel) and the nearest target value in the source aggregate.\n",
    "\n",
    "A powerful feature of `proximity` is that you can target specific values in the aggregate for distance calculation, while others are ignored.  Play with the `target_values` parameter below and see the difference of using `target_values=[1,2,3]` vs. `target_values=[2]` vs. `target_values=[3]`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "   ##### Load data and create `ds.Canvas`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.spatial import proximity\n",
    "from datashader.transfer_functions import dynspread\n",
    "from datashader.transfer_functions import set_background\n",
    "\n",
    "df = pd.DataFrame({\n",
    "   'x': [-13, -11, -5,4, 9, 11, 18, 6],\n",
    "   'y': [-13, -5, 0, 10, 7, 2, 5, -5]\n",
    "})\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H,\n",
    "                x_range=(-20, 20), y_range=(-20, 20))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "   ##### Create Proximity Aggregate\n",
    "   \n",
    "   - Use `Canvas.points` to create an `xarray.DataArray`\n",
    "   - Calculate proximity to nearest non-nan / non-zero elements using `datashader.spatial.proximity`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points_agg = cvs.points(df, x='x', y='y')\n",
    "points_shaded = dynspread(shade(points_agg, cmap=['salmon',  'salmon']),\n",
    "                          threshold=1,\n",
    "                          max_px=5)\n",
    "set_background(points_shaded, 'black')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create proximity grid for all non-zero values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proximity_agg = proximity(points_agg)\n",
    "\n",
    "stack(shade(proximity_agg, cmap=['darkturquoise', 'black'], how='linear'),\n",
    "      points_shaded)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_agg = cvs.line(df, x='x', y='y')\n",
    "line_shaded = dynspread(shade(line_agg, cmap=['salmon',  'salmon']),\n",
    "                          threshold=1,\n",
    "                          max_px=2)\n",
    "set_background(line_shaded, 'black')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "line_proximity = proximity(line_agg)\n",
    "stack(shade(line_proximity, cmap=['darkturquoise', 'black'], how='linear'),\n",
    "      line_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Transform Proximity DataArray\n",
    "Like the other Datashader spatial tools, the result of `proximity` is an `xarray.DataArray` with a large API of potential transformations.\n",
    "\n",
    "Below is an example of using `DataArray.where()` to apply a minimum distance and maximum distance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "where_clause = (line_proximity > 1) & (line_proximity < 1.1)\n",
    "proximity_shaded = shade(line_proximity.where(where_clause), cmap=['darkturquoise', 'darkturquoise'])\n",
    "proximity_shaded = set_background(proximity_shaded, 'black')\n",
    "stack(proximity_shaded, line_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Viewshed\n",
    "\n",
    "The `datashader.spatial.viewshed` function operates on a given aggregate to calculate the viewshed (the visible cells in the raster) for the given viewpoint (observer) location.  \n",
    "\n",
    "The visibility model is the following: Two cells are visible to each other if the line of sight that connects their centers does not intersect the terrain. If the line of sight does not pass through the cell center, elevation is determined using bilinear interpolation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Simple Viewshed Example\n",
    "\n",
    "- The example below creates a datashader aggregate from a 2d normal distribution.\n",
    "- To calculate the viewshed, we need an observer location.\n",
    "- This location is indicated by the orange point in the upper-left of the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.spatial import proximity\n",
    "from datashader.spatial import viewshed\n",
    "\n",
    "from datashader.transfer_functions import dynspread\n",
    "from datashader.transfer_functions import set_background\n",
    "\n",
    "OBSERVER_X = -12.5\n",
    "OBSERVER_Y = 10\n",
    "\n",
    "H = 400\n",
    "W = 400\n",
    "\n",
    "canvas = ds.Canvas(plot_width=W, plot_height=H,\n",
    "                   x_range=(-20, 20), y_range=(-20, 20))\n",
    "\n",
    "normal_df = pd.DataFrame({\n",
    "   'x': np.random.normal(.5, 1, 10000000),\n",
    "   'y': np.random.normal(.5, 1, 10000000)\n",
    "})\n",
    "normal_agg = canvas.points(normal_df, 'x', 'y')\n",
    "normal_agg.values = normal_agg.values.astype(\"float64\")\n",
    "normal_shaded = shade(normal_agg)\n",
    "\n",
    "observer_df = pd.DataFrame({'x': [OBSERVER_X], 'y': [OBSERVER_Y]})\n",
    "observer_agg = canvas.points(observer_df, 'x', 'y')\n",
    "observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),\n",
    "                            threshold=1, max_px=4)\n",
    "\n",
    "normal_illuminated = dsgeo.hillshade(normal_agg)\n",
    "normal_illuminated_shaded = shade(normal_illuminated, cmap=['black', 'white'], \n",
    "                                  alpha=128, how='linear')\n",
    "\n",
    "stack(normal_illuminated_shaded, observer_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Calculate viewshed using the observer location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Will take some time to run...\n",
    "%time view = viewshed(normal_agg, x=OBSERVER_X, y=OBSERVER_Y)\n",
    "\n",
    "view_shaded = shade(view, cmap=['white', 'red'], alpha=128, how='linear')\n",
    "\n",
    "stack(normal_illuminated_shaded, observer_shaded, view_shaded)                         "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Viewshed on Terrain\n",
    "\n",
    "- Let's take the example above and apply it to our terrain aggregate.\n",
    "- Notice the use of the `observer_elev` argument, which is the height of the observer above the terrain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.spatial import viewshed\n",
    "\n",
    "W = 600\n",
    "H = 400\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20e6, 20e6), y_range=(-20e6, 20e6))\n",
    "terrain = dsgeo.generate_terrain(cvs)\n",
    "terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')\n",
    "\n",
    "illuminated = dsgeo.hillshade(terrain)\n",
    "\n",
    "OBSERVER_X = 0.0\n",
    "OBSERVER_Y = 0.0\n",
    "\n",
    "observer_df = pd.DataFrame({'x': [OBSERVER_X],'y': [OBSERVER_Y]})\n",
    "observer_agg = cvs.points(observer_df, 'x', 'y')\n",
    "observer_shaded = dynspread(shade(observer_agg, cmap=['orange']),\n",
    "                            threshold=1, max_px=4)\n",
    "\n",
    "stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "      terrain_shaded,\n",
    "      observer_shaded)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time view = viewshed(terrain, x=OBSERVER_X, y=OBSERVER_Y, observer_elev=100)\n",
    "\n",
    "view_shaded = shade(view, cmap='fuchsia', how='linear')\n",
    "stack(shade(illuminated, cmap=['black', 'white'], alpha=128, how='linear'),\n",
    "      terrain_shaded,\n",
    "      view_shaded,\n",
    "      observer_shaded)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fuchsia areas are those visible to an observer of the given height at the indicated orange location."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zonal Statistics\n",
    "\n",
    "Zonal statistics allows for calculating summary statistics for specific areas or *zones* within a datashader aggregate. Zones are defined by creating an integer aggregate where the cell values are zone_ids.  The output of zonal statistics is a Pandas dataframe containing summary statistics for each zone based on a value raster.\n",
    "\n",
    "Imagine the following scenario:\n",
    "- You are a hiker on a six-day-trip.\n",
    "- The path for each day is defined by a line segement.\n",
    "- You wish to calculate the max and min elevations for each hiking segment as a Pandas dataframe based on an elevation dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.colors import Set1\n",
    "\n",
    "W = 800\n",
    "H = 600\n",
    "\n",
    "cvs = ds.Canvas(plot_width=W, plot_height=H, x_range=(-20, 20), y_range=(-20, 20))\n",
    "\n",
    "terrain = dsgeo.generate_terrain(cvs)\n",
    "terrain_shaded = shade(terrain, cmap=Elevation, alpha=128, how='linear')\n",
    "\n",
    "illuminated = dsgeo.hillshade(terrain)\n",
    "illuminated_shaded = shade(illuminated, cmap=['gray', 'white'], alpha=255, how='linear')\n",
    "\n",
    "zone_df = pd.DataFrame({\n",
    "   'x': [-11, -5, 4, 12, 14, 18, 19],\n",
    "   'y': [-5, 4, 10, 13, 13, 13, 10],\n",
    "   'trail_segement_id': [11, 12, 13, 14, 15, 16, 17]\n",
    "})\n",
    "\n",
    "zones_agg = cvs.line(zone_df, 'x', 'y', ds.sum('trail_segement_id'))\n",
    "zones_shaded = dynspread(shade(zones_agg, cmap=Set1), max_px=5)\n",
    "\n",
    "stack(illuminated_shaded, terrain_shaded, zones_shaded)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datashader.spatial import zonal_stats\n",
    "\n",
    "zonal_stats(zones_agg, terrain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Calculate custom stats for each zone"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = dict(elevation_change=lambda zone: zone.max() - zone.min(),\n",
    "             elevation_min=np.min,\n",
    "             elevation_max=np.max)\n",
    "\n",
    "zonal_stats(zones_agg, terrain, stats)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the zones are defined by line segments, but they can be any spatial pattern, and in particular can be any region computable as a Datashader aggregate.\n",
    "\n",
    "\n",
    "### References\n",
    "- Burrough, P. A., and McDonell, R. A., 1998. Principles of Geographical Information Systems (Oxford University Press, New York), p. 406.\n",
    "- Making Maps with Noise Functions: https://www.redblobgames.com/maps/terrain-from-noise/\n",
    "- How Aspect Works: http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-aspect-works.htm#ESRI_SECTION1_4198691F8852475A9F4BC71246579FAA"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
